Research Article 



Statistics 

in Medicine 



Received 24 May 201 1, Accepted 17 April 2012 Published online 27 June 2012 in Wiley Online Library 

(wileyonlinelibrary.com) DOI: 10.1002/sim.5436 

Misspecification of Cox regression 
models with composite endpoints 

Longyang Wu and Richard J. Cook*^ 

Researchers routinely adopt composite endpoints in multicenter randomized trials designed to evaluate the 
effect of experimental interventions in cardiovascular disease, diabetes, and cancer. Despite their widespread 
use, relatively little attention has been paid to the statistical properties of estimators of treatment effect based 
on composite endpoints. We consider this here in the context of multivariate models for time to event data in 
which copula functions link marginal distributions with a proportional hazards structure. We then examine the 
asymptotic and empirical properties of the estimator of treatment effect arising from a Cox regression model for 
the time to the first event. We point out that even when the treatment effect is the same for the component events, 
the limiting value of the estimator based on the composite endpoint is usually inconsistent for this common value. 
We find that in this context the limiting value is determined by the degree of association between the events, the 
stochastic ordering of events, and the censoring distribution. Within the framework adopted, marginal methods 
for the analysis of multivariate failure time data yield consistent estimators of treatment effect and are therefore 
preferred. We illustrate the methods by application to a recent asthma study. Copyright © 2012 John Wiley & 
Sons, Ltd. 
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1. Introduction 

Many diseases put individuals at elevated risk for a multitude of adverse clinical events, and researchers 
routinely design randomized clinical trials to evaluate the effectiveness of experimental interventions for 
the prevention of these events. Trials in cardiology, for example, record times of events such as non-fatal 
myocardial infaction, non-fatal cardiac arrest, and cardiovascular death [1]. In cerebrovascular disease, 
patients with carotid stenosis can be treated with medical therapy or surgery, and trials evaluating their 
relative effectiveness may record endpoints such as strokes ipsilateral to the surgical site, contralateral 
strokes, and death [2]. In oncology, researchers often design trials to study treatment effects on disease 
progression and death [3], but palliative trials of patients with skeletal metastases may be directed at 
preventing skeletal complications including vertebral and non-vertebral fractures, bone pain, and the 
need for surgery to repair bone [4]. In these and many other settings, although interest lies in preventing 
each of the respective events, it is generally infeasible to conduct studies to answer questions about each 
component. 

When one type of event is of greater clinical importance than others, it can be chosen as the basis 
of the primary treatment comparison, and effects on other types of events can then be assessed through 
secondary analyses. When two or more events are of comparable importance, co-primary endpoints can 
be specified, but tests of hypotheses must typically control the experimental type I error rate through 
multiple comparison procedures [5-7] ; these make decision analyses more complex. A seemingly sim- 
ple alternative strategy is to adopt a so-called composite event [8, 9] that is said to have occurred if any 
one of a set of component events occurs. The time of the composite event is therefore the minimum of 
the times of all component events. 
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There are several additional reasons investigators may consider the use of composite endpoints 
in clinical trials. In studies involving a time-to-event analysis, the use of a composite endpoint 
will mean that more events will be observed than would be observed for any particular compo- 
nent. If the same clinically important effect is specified for the composite endpoint and one of 
its components, this increased event rate will translate into greater power for tests of treatment 
effects; at the design stage this translates to a reduction in the required number of subjects or 
duration of follow-up [9-11]. Composite endpoints are routinely adopted through the introduction 
of one or more less serious events, however, which presumably warrants revising the clinically 
important effect of interest. Moreover, we show later that with models featuring a high degree of 
structure, model assumptions may not even be compatible for the composite endpoint and one of 
its components. 

In time-to-event analyses, interest may lie in the effect of an experimental treatment versus standard 
care on the risk of a non-fatal event. This is a common framework in trials of patients with advanced 
diseases where interest lies in improving quality of life through the prevention of complications. In such 
settings, individuals are at considerable risk of death and a competing risks problem arises. Investigators 
often deal with this by adopting a composite endpoint based on the time to the minimum of the non-fatal 
event of interest and death [12, 13]. This strategy leads to an 'event-free survival' analysis that is par- 
ticularly common in cancer where progression-free survival is routinely adopted as a primary endpoint 
[14]. In palliative trials, however, a treatment may not be expected to have an effect of survival, and 
if a non-negligible proportion of individuals die before experiencing the clinical event of interest, this 
analysis can lead to a serious underestimation of the effect of the treatment [10, 15]. 

Recommendations are available in the literature on how to design trials, analyze resultant data, and 
report findings when composite endpoints are to be used [10-12, 16]. The main recommendations 
include that (i) individual components should have similar frequency of occurrence, (ii) the treatment 
should have a similar effect on all components, (iii) individual components should have similar impor- 
tance to patients, (iv) data from all components should be collected until the end of trial, and (v) 
individual components should be analyzed and reported separately as secondary endpoints. The first 
three recommendations have face validity and seem geared towards helping ensure that conclusions 
regarding treatment effects on the composite endpoint have some relation to treatment effects on the 
component endpoints, thus helping in the interpretation of results. The collection of data on the occur- 
rence of the component endpoints until the end of the trial facilitates separate assessment of treatment 
effects on each of the component endpoints. This means the consistency of findings across components 
can be empirically assessed. 

The aforementioned issues have been actively debated in the medical literature [11, 16-19], but 
there has been relatively little formal statistical investigation of these points. In this paper, we dis- 
cuss statistical considerations related to composite endpoint analyses and use the recommendations 
to guide the investigation. Because the Cox regression model is routinely adopted for the analy- 
sis of composite endpoints in clinical trials [12], we consider it here and point out important issues 
regarding model specification and interpretation. We formulate multivariate failure time models with 
proportional hazards for the marginal distributions that may be used to reflect the settings where 
composite endpoints are most reasonable according to the current guidelines. We study the asymp- 
totic and empirical properties of estimators arising from a composite endpoint analysis. We also 
explore the utility of marginal methods based on multivariate failure time data [20]. We argue that 
the belief that composite endpoints provide an overall measure of the effect of treatment is overly 
simplistic, and a thoughtful interpretation of intervention effects based on composite endpoints alone 
is difficult. Their use as a primary basis for treatment comparison in clinical trials therefore warrants 
careful consideration. 

The remainder of this paper is organized as follows. In Section 2, we construct bivariate failure 
time distributions for which the marginal distributions have proportional hazards between two treat- 
ment groups. We then derive the distribution for the time to the first event and show that it does not 
typically feature proportional hazards across the two treatment groups. We use large sample theory for 
misspecified models to derive the limiting value of the log hazard ratio from a naive Cox model, and 
empirical studies demonstrate finite sample properties which are in close alignment with the theory. 
An alternative approach to synthesizing data over component events is to conduct a global analysis on 
the basis of the marginal methods of Wei et al. [20]; we explore this in Section 3. An application to a 
recently completed asthma management study illustrates the various methods in Section 4, and we make 
the concluding remarks in Section 5. 
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2. Multivariate failure time distributions via copula functions 

2.1. Construction of joint distributions based on copula functions 

If (Ui, U 2 )' is a bivariate random variable with standard uniform margins on [0, 1], a two-dimensional 
copula function can be defined as 

C(u u u 2 ) = P(U l ^ui,U 2 ^u 2 ) , (1) 

[21]. If there exists a convex decreasing function T-L(u; 9) such that T~L : (0, 1] — ► [0, 00) and T-L{\ ; 9) — 0, 
and if the copula function can be written as 

C(ui,« 2 ; 0) = n^CHiu^ 9 ) + H(u 2 ; 9); 9) , 

then copula belongs to the Archimedean family; the univariate function T-l{u\9) is called the generator 
for the copula [22]. Suppose (Un, Ui 2 )' and (Uj\, Uj 2 )' are two random variables drawn from the joint 
distribution (1). A common measure of the association between U\ and U 2 is Kendall's r, defined as 

re = P {(U n - U n )(U i2 - U j2 ) >0;9}-P {(U n - U n )(U i2 - U j2 ) < 0; 9} . 

where we write tq to make the relation between 9 and r explicit. 
For Archimedean copulas, this can be written as 

f 1 U{u;9) 
r e = l+4/ — — — -du. 
Jo U'{u;9) 

Copula functions have received considerable attention in the statistical literature in the past few years 
because they offer a convenient and attractive way of linking two marginal distributions to create a joint 
survival function [23]. Suppose T\ and T 2 are a pair of non-negative random variables with respec- 
tive survivor functions ^(t^z^i) and J- 2 (t 2 \z; ct 2 ) given a covariate z. If we let U\ — F\{T\\z;ot\) 
and U 2 — J- 2 (T 2 \z; u 2 ) where indexes the marginal distribution for T^z, then Uk ~ UNIF(0, 1), 
k — 1,2. We can define the bivariate 'survival' distribution function of (U\, U 2 ) through a copula as in 
(1) and obtain a joint survivor function for (T\, T 2 )' given Z as 

Ti 2 (h,t 2 \z; Q) = PiTr Zh,T 2 ^ t 2 \z; tt) = CgiTiihlz;^), F 2 (t 2 \z;ce 2 ); 9) , (2) 

where £2 — (a'. 9)' with a — (a^a^)' ■ Because Kendall's x is invariant to monotonic increasing 
or decreasing transformations [21], it can be interpreted as a measure of association of the trans- 
formed variables (T\, T 2 )' given Z. The use of a copula function to define the joint distribution of 
{T\ , T 2 )\z is particularly appealing because one can specify the marginal distributions to have a propor- 
tional hazards form; this is not typically possible for joint distributions induced by random effects or 
intensity-based analyses. 

If a composite endpoint analysis is planned, it would be based on modeling the random variable 
T — min(7i , T 2 ), which has survival, density, and hazard function conditional on Z, given by 

P(T Zt\z) = T(t\z\ £2) = Fi 2 (t,t\z\ £2) , (3) 

f(t\z) — — dJ^lz; £2)/d/, and X(t\z; £2) = — dlog ,F(/|jc; £2)/dt, respectively. Suppose Z is a binary 
indicator where Z = 1 for individuals in a treatment group and Z —0 otherwise. A key point is that the 
hazard ratio X(t\z — 1; Q)/h(t\z = 0; £2) is not, in general, independent of time. As a result, even if the 
marginal distributions feature proportional hazards, the model for the composite endpoint will typically 
not. We study this point further in the next four settings for three different Archimedean copulas and the 
case of independent components. 
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2.7.7. Composite endpoint analysis based on a Clayton copula. The Clayton copula [24] is a member of 
the Archimedean family with generator 'H (u; 9) — u~ s — l,7^ _1 (i;; 9) — (v + \)~ l l e and copula function 



C(ui,u 2 ; 0) = (u~ 9 + u~ 9 - i\ 



-lie 



(4) 



with 9^—1. Kendall's x is then given by tq — 9/(9 + 2), which can be seen to vary over [—1, 1]. 

Consider the joint distribution of (T\,T2)\Z in which the marginal distribution for T k \Z, 1,2 
features proportional hazards; so X k (t\z) — X k o(t) expifikZ) with A k (t\z) — A k0 (t)exp(fi k z) where 
A k0 (t) — L X k0 (s)ds, k — 1,2. If the joint survivor function Tn{t\, t%\z\ ^) is determined by the 
Clayton copula through (2), by (3) the survivor function of the failure time T — min(ri , T 2 ) given z is 



F{t\z\ O) = [exp(9A 10 (t)e PlZ ) + exp(0A 2 o(Oe fcz ) - l] 



Hence, the hazard ratio for the treatment versus control groups for the composite endpoint is 



X(t\z. = 1;S2) 


ELi x ko(t) exp (p k + 9A k0 (t)e^) 


1 


ELi exp(0A fco (Oe fe ) - 1 


X(t\z. = 0; Q) 


J2l=iho(t)exp(9A ko (t)) 


1 


ELi ex P( 0A feo(O)-i] 



(5) 



(6) 



which is not invariant with respect to time in general. 

To gain some insight into this function, suppose the marginal distributions are exponential with 
common baseline hazards of A 10 (0 = A 2 o(f) — X — log 10 so that the probability of a type k event 
occurring before t — 1 is 0.90 for a control subject (i.e., P(T k < 1|Z = 0) = 0.90). Further suppose 
that a common hazard ratio of 0.50 holds for the two margins (i.e., exp(/3i) = exp(^2) = 0.50). This 
setting is consistent with the recommendations that the component events occur with comparable fre- 
quency because P{T\ < Ti\Z) — 0.5, and have comparable treatment effects — fc)- Figure 1(a) 
contains a plot of the hazard ratio (6) over the time interval [0, 1] for models with mild (tq — 0.2), mod- 
erate (tq — 0.40), and strong (tq — 0.60) association. As can be seen, even when the treatment effects 
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Figure 1. Plots of the hazard ratio over the time interval [0, 1] for the composite endpoint model implied 
by the Clayton copula (panel (a)) and Frank copula (panel (b)) with marginal exponential distributions with 
Xi = X 2 = log 10 and exp(0i) = exp(/J 2 ) = exp(/J) = 0.50 and mild (xg = 0.20), moderate (re = 0.40), and 

strong (tq = 0.60) associations. 
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are the same for the two component endpoints, there can be non-negligible variation in the hazard ratio 
over time, and within this family of models, the nature of this variation depends on the strength of the 
association between the two failure times. 



2.7.2. Composite enclpoint analysis based on a Frank copula. The generator for the Frank copula [25] 
is H(u; 9) — — log((exp(— 9t) — l)/(exp(#) — 1)), and the resulting copula function is 



C(wi, u 2 ; 9) — 



los 



(e 



-9ui 



iu 2 



1) 



- 1 



where 9 e 5R; Kendall's t is then t e = 1 - A9~ l + 49~ 2 f Q r/(exp(r) - l)df. If we adopt the same 
marginal distributions as before, the survivor function for the composite endpoint is 



nt\z) 



log 



(exp(-6»e- Al(0 ^ lZ )- 1) (exp(- 



-A 2 {t)e e 2' 



)-l) 



1 



but because A(f |z; £2) = — dlog J 7 (t)/dt, the hazard ratio X(t\z — 1; £2)/X(t\z — 0; Q) has a compli- 
cated form. Figure 1(b) contains a plot of this hazard ratio over [0, 1], and as in the case of the Clayton 
copula, there is considerable variation in this ratio over time. 



2.1.3. Composite endpoint analysis based on a Gumbel-Hougaard copula. The generator for the 
Gumbel-Hougaard [26] copula is 'H(u\ 9) — (—log?) 6 * giving 



C( Ml ,M 2 ;0) = exp(-((-log Ul ) e + (-logM 2 )y , 

fox9 > 1; Kendall's x is given by tq — (9— \)/9. The corresponding survivor function for the composite 
endpoint is 



F(t\z) = vq?\ 
and if fix — fi 2 — fi, the hazard is 

k(t\z) = exp(Pz) 



(A 1 (0^ lZ ) e + (A 2 (0e foz ) ( 



(A!(f) e + A 2 (O e ) 6 



Ai(0Ai(0 6 



A 2 (/)A 2 (0 6 



Interestingly, the hazard ratio in this case is exp(/3), which means that the proportional hazards model for 
the composite endpoint is compatible with a proportional hazards model for the margins. If the hazard 
ratio is in fact common for the component endpoints, then a consistent estimator will be obtained for this 
common effect on the basis of a Cox model for the composite endpoint. 



2.1.4. Composite endpoint analysis with independent components. Here, we consider the setting 
where the component failure times are independent; a special case of tq — 0 for the joint models in 
Sections 2.1.1-2.1.3. In this case, the hazard ratio for the composite endpoint analysis reduces to 

Ht\z = l; a) _ A 10 (?) exp(j8i) + A 20 (Q exp(^ 2 ) 
X(t\z = 0;a) Aio(0 + A 20 (0 

With nonhomogeneous hazards, it is apparent that the composite endpoint analysis is only compatible 
with a proportional hazards assumption if either (A.l) j8i = /J 2 = fi or (A.2) Ai 0 (0 = A 20 (0- If 
jf3i = ($2 — fi, then a consistent estimate of this common effect is obtained in a composite endpoint anal- 
ysis. If yf3 1 ^ /3 2 but the hazard functions are identical, the multiplicative effect is (exp(/3i) + exp(/3 2 ))/2. 
If assumptions A.l and A.2 do not hold, then the ratio is a complicated time varying function of the 
baseline hazards and respective treatment effects. 
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2.2. Misspecification of the Cox model with composite endpoints 

The previous section demonstrated that the composite endpoint analysis is typically based on a misspec- 
ified Cox regression model if the marginal distributions satisfy the proportional hazards assumption. In 
this section, we investigate the frequency properties of estimators from a composite endpoint analysis 
when the component endpoints are associated through a copula function. 

Let Ti — min(7h , T^) denote the time of the composite endpoint for individual i in a sample of size 
m. Let {Ni(s),s < 0} denote the counting process for subject i, which indicates the occurrence of the 
composite endpoint, so that dNi(s) = 1 if 7* = s and is zero otherwise. Suppose that it is planned 
to follow all subjects over the interval (0, C'] but that subjects may be lost to follow-up or withdraw 
from the study prematurely. Let Wi represent the withdrawal time for subject i and Q = min(M^-, C*) 
denote their right censoring time. Let Yj(s) — I(s ^ Ti) indicate whether subject i is at risk of the 
composite endpoint at time s, Y^(s) — I(s ^ C,) indicate whether they are under observation at time 
s, and Yi(s) — Y?(s)Yi(s) indicate whether they are event free and under observation. The observ- 
able counting process for the response is then based on dNj(s) — Yi(s)dNi(s) for subject i. The 
data for a sample of size m then consist of {(Yj(s),dNi(s)) ,0 < s Zj,i — 1, . . . ,m}, which if we let 
Y(s) = (fiCi), . . . , Y„(s)y, d_N(s) = (dNi(s), . . . ,dN m (s))' and Z = (Z u . . . , Z m )', we may write 
more compactly as {(y(^), dN(s)) , 0 < s, Z}. 

The Cox model is widely used in the analysis of composite endpoints [27] to estimate the relative 
hazard where we assume the hazard function for Ti \n to have the form 

f(t\Zi) - fo{t) exp(azi), (7) 

where ^o(0 is a non-negative baseline hazard function corresponding to the control group and z.i is the 
treatment covariate for individual i,i — 1, . . . ,m. The treatment effect a can be estimated using the 
maximum partial likelihood [28] by solving 

A f°° ( S^(a t)\ 
U(a) = J2J 0 fi « ( z > '■- I d ^ (0 (8) 



i = i 



where S (k) (a, t) = >ZT=i Yi(t)zf exp{az;}, k = 0, 1. 

If |y/(^),0 < s\ is independent of {Ni(s),0 < s} given Z, and if (7) is correctly specified, then 

(8) has expectation zero and the solution a is consistent for the true value, a. In the independence 
case, this true value is /3 if the treatment effect is common (i.e., under A.l j3 — fi\ — fa) or 
a — log(exp(/3i) + exp(/32))/2 if the baseline hazard functions are the same (i.e, under A. 2). More 
generally, however, a is consistent for a*, the solution to expected score function U(a) — E(U(a)) 
given by 

U{a) = C \ E (^ f<(0 ^ <(0 )" g(5w[alo) g (E f 'W<WK0j|, (9) 

where the expectation E is with respect to the true model for dN(s)) , 0 < s, Z} [29-31]. By 

using the true model based on (3) and assuming independent censoring for the withdrawal time Wi with 
survival distribution Q(w\z) — Q(w), these expectations can be obtained as follows: 

E(S (1) (a, t j) = m[G(t)F(t\Z; Q)] exp(a)P(Z = 1) 

£(5 (0) (a, 0) = m[g(t)T(t\Z = 1; Q)] exp(a)P(Z = 1) + m\Q(t)T(t\Z = 0; Q)]P(Z = 0) . 



Likewise, 



E (J2™ =1 Y i(t)dNi(t)) = mQ{t) Y^ r=0 ?(t\Z = r; Q)X(t\Z = r)P(Z = r), 
E (J2™ =1 ZiYi(t)dNi(t)) = mQ(t)\T(t\Z = 1; Q)X(t\Z = l)P(Z = 1)] . 
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To illustrate the bias resulting from a composite endpoint analysis, consider a randomized 
clinical trial in which subjects are to be followed over the interval (0, C *] where C* = 1. Let Z = 1 
for treated subjects and Z = 0 for control subjects and suppose P(Z = 1) = 1 — P(Z — 0) = 0.5. 
We set j6i = /J 2 — P — log 0.80 to consider the case compatible with the current recommenda- 
tions on the use of composite endpoints. We set Ai and A 2 so that (i) P(Ti < T 2 \Z — 0) = p\ 
equals a desired probability that the type 1 event occurs before the type 2 event among control sub- 
jects and that (ii) P(C^ < T) — it a satisfies the administrative censoring rate for the composite 
endpoint among all subjects, where jia — 0.20. Finally, suppose subjects may withdraw from the 
study early, and let W have an exponential distribution with rate p such that P(C < T) — n, where 
P(C < T) = E Z [P(W < T < C f |Z) + J P(C t < T\Z)] and n is the overall censoring rate set to 
it = 0.20, 0.40, 0.60, and 0.80. 

Figure 2 shows the limiting percent relative bias (100 (a* — fi) //?) of the treatment coefficient from 
a composite endpoint analysis when the data are generated by a Clayton copula with mild (r = 0.20) 
and moderate (t = 0.40) association. We plotted this relative bias against P(T\ < Ti\Z = 0) = pi, 
and interestingly, the bias is greatest when p\ — 0.50 but decreases as this probability approaches zero 
or one. In either of the extreme cases {p\ — 0 or p\ — 1), the composite endpoint coincides with the 
occurrence of a single endpoint, and a consistent estimate of the common treatment effect is obtained. 
Note that the bias (a* — /J) is positive, and hence, the limiting value of the treatment effect is more 
conservative than the true common value for each of the components. This means that the estimated 
value would, on average, under-represent the magnitude of the treatment effect on either component, a 
conclusion in line with the findings of [10, 15]. Moreover, we note that the common event rate and the 
common treatment effect are precisely the setting where composite endpoints are recommended for use 
[10-12, 16]. The plots also reveal the sensitivity of the limiting value to the degree of random censoring; 
the higher the censoring rate, the smaller the asymptotic bias. This highlights an important point that the 
limiting value of an estimator from a misspecified failure time model is highly sensitive to the censoring 
distribution even under independent censoring. By comparing the left and right panels in Figure 2, it 



'J i 




Figure 2. Asymptotic percent relative bias (100(a* — ft)/ ft) of Cox regression coefficient of treatment effect 
from composite endpoint analysis when bivariate failure times are generated by a Clayton copula; exponential 
margins, 20% administrative censoring (tta = 0.20), 50:50 randomization, exp(/?i) = exp(/32) = 0.80, and four 
different degrees of additional random censoring (none, 20%, 40%, and 60%). 
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is also apparent that the asymptotic bias is dependent on the degree of association between T\ and T 2 ; 
the greater the association, the greater the asymptotic bias. This makes sense because when the event 
times are independent, consistent estimates should be obtained because assumptions A.l and A. 2 of 
Section 2.1.4 are satisfied. 

Although of secondary interest, one can also show that ifo(t), 0 < t < C^, is consistent for 

= T!r= 0 Q(t\Z - r)T(t\Z = r)X(t\Z = r)P(Z = r) 
° J2 1 r=oS(t\Z = r)T(t\Z = r)exp(a*I(r = l))P(Z = r) 

which when P(Z — 1) = 0.5 and the censoring distribution is the same in the two groups reduces to 

f*(t) = \T(t\Z = \)\{t\Z = 1) + T(t\Z = 0)X(t\Z = 0)]/[F(t\Z = 1) exp(o!*) + F(t\Z = 0)] 



2.3. Simulation studies involving composite endpoints 

2.3.1. Simulation design. Here, we simulate data from (2) to examine the empirical performance of 
estimators for finite samples. We assume that given Z, has an exponential distribution with hazard 
Afc exp(/3fcZ), k — 1,2, and model the association between T\ and T 2 through a Clayton copula. We 
let T — min(ri, T2) denote the time of the composite endpoint as before. We suppose interest lies in 
following subjects over (0, 1]. As in the previous section, we determined the parameters Ai and A2 to 
satisfy the constraints P(T\ < T 2 \Z = 0) = p\, where p\ — 0.25, P(C^ < T) — jta, and we set 
the administrative censoring rate to it a — 0.20. We also incorporated random loss to follow-up with an 
exponential withdrawal time giving a net censoring rate of tc — 0.20, 0.40, 0.60, and 0.80 subject to the 
constraint it a ^ ft- 

For each parameter configuration, we derived the sample size for the composite endpoint analysis 
to achieve a prespecified power under the assumption that the Cox model in (7) holds. Therneau and 
Grambsch (2000) show that the required number of events is D — A{z.\- Yl + Zi_y 7 ) 2 /(a*) 2 , where z q 
is the qth quantile of the standard Normal distribution, yi is the type I error for a one-sided test, 1 — y 2 is 
the power, and a* is the limiting value of treatment effect estimate obtained from (7). We focus on two- 
sided tests at the 5% significance level (yi = 0.05) and sample sizes to achieve 80% power (y 2 = 0.20). 
We calculated the required number of subjects as m — D/P(T < C). In all simulation studies, we con- 
sidered both equal treatment effects — j3 2 — — —.223) and unequal treatment effects (fii — —.223 
and fl 2 — 0)- For eacn parameter configuration, we generated 2000 replicates. We report the mean of the 
a estimates, the empirical standard error (ESE), the average model-based standard error (ASE\), and the 
average robust standard error (ASE 2 ). We also reported the empirical coverage probability (ECP*%) of 
nominal 95% CIs for a* based on robust standard errors and the empirical coverage probability of these 
intervals for /3j (ECP%). The last column contains the empirical power (EP%) of a Wald test of the null 
hypothesis of no treatment effect. 



2.3.2. Composite endpoints with dependent components. Table I contains the simulation results with 
dependent component times given by x — 0.40. The results for equal treatment effects are given in the 
top half of the table that we comment on first. The fourth column contains a*, the limiting value of the 
estimator from the misspecified Cox model in (7). The fact that these values are all smaller in absolute 
value than the true common effects reveals the conservative nature of this limiting value, as already dis- 
cussed in relation to Figure 2; the dependence of the limiting value on the degree of censoring is also 
apparent. This limiting value was used to derive the sample size (m) in the third column. The average 
estimator from the fitted Cox models reported in the fifth column closely approximates the limiting value. 
There is also close agreement between the empirical, average model-based, and average robust standard 
errors. The empirical coverage probabilities of the robust 95% CIs are very close to the nominal levels, 
and the empirical power is in good agreement with the nominal power of 80%. It is worth emphasizing 
that the empirical coverage probability is computed for the parameter a*, not the common /3; for this 
latter parameter, the coverage rates are considerably lower. 

In the bottom half of Table I, we reported the results for the case fi\ ^ fi 2 , where a* is consider- 
ably smaller than This smaller limiting value leads to considerably larger sample sizes to achieve 
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Table I. Frequency properties of estimators of treatment effect based on 
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7Ta = P(C' < T) is the administrative censoring rate, it = P(C^ < T) is the net censoring rate, ESE is the empirical 
standard error, ASEi is the average model-based standard error, ASE2 is the average robust standard error, ECP*% is 
the empirical coverage probability for a* of nominal 95% CIs using the robust standard error, ECP% is the empirical 
coverage probability for fi\ of nominal 95% CIs using the robust standard error, and EP% is the empirical power of a 
Wald test of Hq :a = 0 based on the robust standard error. 



the desired power. Again, however, we see close agreement between the average estimate and the lim- 
iting value, and very close agreement between the average model-based and average robust standard 
errors. The empirical coverage probability (for a*) is also consistent with the nominal level, as is the 
empirical power. 



2.3.3. Composite endpoints with independent components. Table II presents the simulation results with 
independent components (i.e., x — 0). The results in the top half of Table II reveal that the limiting value 
a* is the same as the common value [3 — 0i = 02 as expected because assumption A.l of Section 2 
is satisfied. Again, the average point estimate is in close agreement with this common value, and the 
three standard errors are in close agreement. When the treatment has an effect on Ti and not T2, a* 
is again considerably smaller than [3\. Note, however, even though this is a misspecified model, the 
limiting value does not depend on the censoring distribution. This much smaller value leads to larger 
sample size requirements than in the top half of the table. Because the first component T\ happens less 
frequently than the second component T2 (i.e., P(T\ < 72 |Z = 0) = 0.25), the limiting value from the 
misspecified Cox model is heavily attenuated in this setting. However, neither administrative nor random 
censoring appears to affect the limiting value of the estimator of treatment effect. 



3. A multivariate semiparametric analysis 

3.1. Limiting values for a Wei-Lin-Weissfeld analysis 

In this section, we investigate the utility of the marginal approach of Wei et al. [20] for handling multi- 
variate failure time data. This approach is based on formulating ordinary Cox models for each component 
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Table II. Frequency properties of estimators of treatment effect based on 
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7Ta = P(C^ < T) is the administrative censoring rate, n = P(C^ < T) is the net censoring rate, ESE is the empirical 
standard error, ASEi is the average model-based standard error, ASE2 is the average robust standard error, ECP*% is 
the empirical coverage probability for a* of nominal 95% CIs using the robust standard error, ECP% is the empirical 
coverage probability for fi\ of nominal 95% CIs using the robust standard error, and EP% is the empirical power of a 
Wald test of Hq : a = 0 based on the robust standard error. 



event to obtain component-specific estimates of treatment effect, and it is therefore compatible with the 
way the joint distributions were constructed using copula functions in Section 2. Estimation proceeds 
under a working independence assumption, as often adopted for analyses based on generalized estimat- 
ing equations. We obtain a robust estimate of the covariance matrix, and then we obtain a global estimate 
of treatment effect by taking a weighted average of all component-specific estimates with weights cho- 
sen to minimize the variance of the global estimator. A key distinction between the global approach of 
Wei et al. [20] and the composite endpoint approach is that the former makes use of all observed events 
whereas the composite endpoint uses only information on the first event. 

In the derivations that follow, the composite endpoint is comprised of K components, but we subse- 
quently focus on the case K — 2. We let dNi k (s) — I{Ti k — s), {Nf k (s), 0 < s} denote the counting 
process for type k events, and {Ni(s) — (Nn(s),Ni2(s), 0 < s} denote the bivariate counting process 
for subject i, i = 1, . . . ,m. Let Y ik (s) = I(s s= T ik ), Y?(s) = I(s s= Q), and Y ik (s) = Y?(s)Y ik (s), 
k = 1, . . . , K, i = 1, . . . ,m. The Cox model for the type k event is 

where X k0 (t) is the baseline hazard function and f} k is the corresponding treatment effect. The k\h 
component- specific score function for /3 k is 

m f°° - ( S^(Bh t)\ 
U k (fa) = Y f ik (t) Zi - k ( a J ] dN,k(f) , (10) 

t \ m - 

where >tf, u) = £ Y ik (t)z r t exp{^}, r = 0, 1. 

i = l 

Under the copula model of Section 2 with marginal distributions featuring proportional hazards, 
the solution to the score equation (10), fi k , is consistent for the true treatment effect /3 k . If we let 
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/} = (jSj, . . .,f} K )' and its estimate j} = (pi, . .. , Pk) T , Wei et al. [20] show that y'm (/?-/?) converges 
in distribution to a multivariate normal distribution with a zero-mean vector and variance-covariance 
matrix £ and provides a consistent sandwich-type estimate for £ (/J). 

The global estimate of treatment effect proposed by Wei et al. [20] is a linear combination of all 
component- specific treatment effect estimates fii , . . . , Pk and can be obtained as 

?=c(M. (id 

where the weight c(/3) = £ J[J'£ (/3) _1 J] -1 is chosen to estimate the weight matrix to minimize 
the variance in the class of all linear estimators; £ (/3) is the estimate for the variance-covariance matrix 
of/3andJ = (l,...,l)'. 

To compare the performances of the global approach and the composite endpoint approach, we obtain 
the limiting value of f> as 

P* = c(fi)'p, (12) 

where c(/3) = T,- 1 (fi)J[J'X- 1 (fi)S\- 1 . We therefore require the limiting value of the robust variance 
£ to obtain the limiting value. The detailed derivations are deferred to the Appendix. 

An alternative asymptotically equivalent approach to estimating the global effect and to deriving the 
limiting value involves specifying a single Cox regression model and fitting it using all events while 
'stratifying' on the event type [32]. Although this has some appeal, we adopt the current framework on 
the basis of synthesizing estimates from separate Cox regression models because it makes explicit the 
fact that the global estimate and associated limiting value may be viewed as a weighted average of the 
component- specific estimates. 



Table III. Empirical properties of the global estimates of treatment effect based on Wei-Lin-Weissfeld 
analysis: data were generated under a Clayton copula with r = 0.40, fii = —0.223. 
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7Ta = P(C^ < T) is the administrative censoring rate, n = P(C < T) is the net censoring rate, ESE is the empirical 
standard error, ASEi is the average model-based standard error, ASE2 is the average robust standard error, ECP*% is 
the empirical coverage probability for f}* of nominal 95% CIs using the robust standard error, ECP% is the empirical 
coverage probability for /Si of nominal 95% CIs using the robust standard error, and EP% is the empirical power of a 
Wald test of Hq : f} = 0 based on the robust standard error. 
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3.2. Comparison of the global approach and the composite endpoint analysis 

Table III reports the results from a global analysis of treatment effect based on the marginal analysis 
proposed by Wei et al. [20]. In this table, the sample sizes were computed on the basis of the formula 
for the composite endpoint analysis using the limiting value of the regression coefficient. As one would 
expect from (10), when the treatment effects are equal, then the marginal analysis yields consistent esti- 
mators for this common effect and the mean estimate across all simulated trials is very close to the 
limiting value. Moreover, the ESE and the average robust standard error were in very close agreement; 
the average model-based standard error is conservative because it is based on the working independence 
assumption being correct. The empirical coverage probabilities (based on the robust standard errors) 
were compatible with the nominal 95% level for fa when fa — fa. When fa ^ fa, the empirical 
coverage for fa was zero, a reflection of the difference between fa and fa. When fa = 0, the limiting 
value fa was quite small, and hence, the sample sizes of the trial were much larger. 

When fa ^ fa, the composite endpoint and global analyses yield estimators that do not coincide with 
fa, fa, or each other. We next compare these limiting values. We consider the case in which two fail- 
ure times are generated by a Clayton copula with exponential margins and a single treatment covariate 
modeled through proportional hazards with fa — log(0.80) and fa — 0. We consider mild and mod- 
erate association between the failure times with r — 0.20 and r — 0.40, respectively. Administrative 
censoring was set to 40% and additional random censoring from an exponential withdrawal time gave 
cases with 60% and 80% as well. The limiting values of the composite endpoint and global analyses 
were plotted against P(T\ < T 2 \Z = 0) = p\ in Figure 3. It is apparent that when p\ approaches 
zero, the limiting value for both methods approaches 0. For the composite endpoint, this makes sense 
because the first event is most likely to be a type 2 event for which there is no treatment benefit. As 
Pi approaches 1, the limiting value for the composite endpoint analysis approaches fa for analogous 
reasons. The limiting value from the global analyses tracks these limiting values quite well, but tend to 
correspond to larger estimates of treatment effect because the limiting value is larger in absolute value. 
Thus, even when the two components have equal frequencies and the proportional hazards assumption 
holds for each component, the global analysis, in the limit, will yield an estimate of treatment effect 
that is greater than that of the composite endpoint analysis. These relationships hold across both lev- 
els of association and over different degrees of censoring. Although we have restricted attention to the 
Clayton copula in these calculations and empirical studies, this investigation could be repeated under 
other copula models, and although the limiting values would differ, qualitatively similar findings would 
be expected. 



4. Application to an asthma management study 

We now apply both the composite endpoint analysis and the global approach to an asthma management 
study [33]. This is a two-phase, multicenter, randomized, parallel group effectiveness trial for comparing 
two treatment strategies for asthma management over a 2-year period. The control strategy is a 'clinical 
strategy', in which the treatment was guided on the basis of patient symptoms and spirometry readings. 
The experimental strategy is a so-called 'sputum strategy' (SS), whereby a cellular analysis of sputum 
samples was used to guide corticosteroid therapy use to keep eosinophils cell counts less than 2%. In 
phase I, a total of 107 patients were identified through the minimum treatment to maintain control. 
The aim of this asthma study was to investigate whether SS is more effective than clinical strategy on 
reducing the number and severity of exacerbations in phase II. 

In our analysis, we focus on two types of exacerbations: mild exacerbations defined as requiring a 
daily maintenance dose of fluticasone of < 250 |ig and severe exacerbations defined here as requiring a 
minimum daily maintenance dose ^ 250 |ig. The composite endpoint is defined as the time to the first 
of the two type of exacerbations. Figure 4 displays the empirical distribution function plots for the two 
component types of exacerbations and for the composite endpoint. It is apparent that the severe exac- 
erbations occur much more frequently than mild exacerbations and thus represent the majority of the 
events contributing to the composite endpoint. 

Table IV presents the results of the proportional hazards regression analysis in which the single binary 
covariate is the treatment indicator taking the value one for patients in the experimental (SS) group 
and zero otherwise. From these results, it is clear that the experimental SS strategy leads to a signifi- 
cantly lower hazard of severe exacerbations with a relative risk reduction of 47% (95% CI: 0.02, 0.71; 
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Figure 3. Plot of limiting values of regression estimates of treatment effect based on a composite endpoint anal- 
ysis and a global Wei et al. [20] analysis with bivariate data generated with a Clayton copula; fi\ = log 0.80, 

p 2 = o. 
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Figure 4. Empirical distribution functions for severe exacerbations, mild exacerbations, and the composite 

endpoint in asthma trial. 
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Table IV. Analysis results of the asthma management study. 



Endpoint/ Analy si s 


RR 


95% CI 


p-value 


p* 


Severe 


0.53 


(0.285, 0.977) 


0.042 


0.22 


Mild 


2.14 


(0.624, 7.310) 


0.227 


0.11 


Composite 


0.665 


(0.388, 1.138) 


0.137 


0.063 


Global (WLW) 


0.702 


(0.405, 1.219) 


0.209 





p =0.042) but has little effect on the occurrence of mild exacerbations (p — 0.227). The result from the 
composite endpoint analysis is not statistically significant with p — 0.137. The Wei-Lin-Weissfeld [20] 
global analysis yields an estimate that is close to that obtained from the composite endpoint analysis, 
but it is apparent when examining the effects on the separate components that a global estimate is not 
an adequate summary of the data. The last column of Table IV gives the /^-values for testing the propor- 
tional hazards assumption using univariate tests based on Schoenfeld residuals [32]. There is insufficient 
evidence to reject the null hypothesis of proportional hazards for each component, and the test yields a 
/7-value just shy of statistical significance for the composite endpoint analysis at 0.063. Thus, although 
we have demonstrated that, in principle, if the proportional hazards assumption holds for the components 
of a composite endpoint, it generally does not hold for composite endpoint itself, the tests do not suggest 
problems with model fit for this particular data. Although the association may be well characterized by 
the Gumbel-Hougaard copula, the power of the tests for departures from proportional hazards may also 
be inadequate. 



5. Discussion 

Composite endpoints are widely adopted in clinical trials, and fitting a Cox proportional hazards model 
is the standard approach to estimating treatment effects on the basis of such endpoints. We have demon- 
strated that even when the treatment effects are the same for component endpoints under marginal Cox 
models, the Cox model for the composite endpoint is typically misspecified because the proportional 
hazards assumption does not in general hold. The estimator of treatment effect under such a misspeci- 
fied Cox model for the composite endpoint has a slightly conservative limiting value, meaning that the 
benefit of treatment was under-estimated in the settings we examined. We found several factors that 
influence the limiting value including the strength of the association between the individual component 
events, stochastic ordering of the individual components, and the degree and nature of the censoring 
process; empirical studies corroborated these findings. Although we have not explored this here, it is 
clear from Section 2 that the specific copula function would also have an important effect. More gener- 
ally, variation in the treatment effect across the individual components makes it even more difficult to 
interpret estimators. 

Composite endpoints are often thought to offer a measure of the 'overall effect' of a treatment [9]. 
In fact, the opposite can be true if treatment effects are in opposing directions for different components, 
and one component tends to occur first. The event tending to occur first will have the greatest influence 
on the estimator of treatment effect based on the composite endpoint, masking the effect on the events 
that tend to occur later. In the case where the events occur with equal frequency and the treatment effect 
is the same for the components, the asymptotic calculations of Section 2 show that the estimate based on 
the composite endpoint suggests a smaller benefit that holds for the components. Although some might 
argue that this is therefore a conservative approach, it is essentially incorrect. The global approach of 
Wei et al. [20] can, however, provide evidence of this adverse effect in the component-wise analysis, and 
this will lead to an attenuation of the global effect in the weighted analysis. 

Another rationale put forward for adopting composite endpoints is to model the event-free survival 
probability. For example, Sheehe [34] proposed that the event-free survival curve can be computed on 
the basis of Cox model estimates of hazard ratios from the composite endpoint containing mortality as 
a component. As we have demonstrated, estimators from Cox regression using composite endpoints can 
be attenuated by including components for which there is no treatment effect. Using estimates from the 
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composite endpoint analysis (event-free survival) may not provide a valid representation of the effect on 
the non-fatal event. 

Two of the guidelines for the use of composite endpoints include the requirement that individual com- 
ponent events should be of roughly equal frequency, and the treatment effects should be comparable 
across all components. Our analytical and empirical investigation shows that these may not be sufficient 
conditions if interest lies in estimating these common effects in the sense that even when these condi- 
tions are satisfied, the association between the two events can lead to substantial bias in estimators based 
on composite endpoints. We support the recommendations that (i) data from all components should 
be followed until the end of the trial and (ii) individual components should be analyzed and reported 
separately. This alternative design and analysis facilitates a global approach [20] based on combining 
estimates from individual components, as well as assessment of whether this is appropriate. In the con- 
text of the copula-based joint model, we found that the global approach, in general, outperforms the 
composite endpoint analysis in terms of the properties of the resulting estimators and power or sample 
size requirements. 

We have formulated a model with proportional hazards for each component event through the use of a 
copula function to reflect an idealized situation in alignment with the recommendations in the literature. 
We restricted attention to the situation with two component endpoints, but three or more components 
are often specified in practice. When multiple components are of interest, copula functions with an 
'exchangeable' association structure can be readily adopted; more baseline marginal hazard functions 
and treatment effects would need to be specified. It is relatively straightforward to extend the derivations 
and empirical studies reported here for this setting but more challenging to cover a meaningful spectrum 
of settings, summarize results, and make recommendations. 

Alternative frameworks could naturally be adopted for specifying models for correlated failure time 
data. One might, for example, consider intensity-based models where the risk of one type of event 
changes with the occurrence of another type of event. This could arise because of a biological mech- 
anism in which the medical risk actually increases or if treating physicians alter the therapy being 
given. This formulation, although natural for characterizing the response process, is not compatible 
with proportional hazards for the marginal models. One might also consider frailty models for address- 
ing the association between event times, but again, the marginal models will not have a proportional 
hazards form. 

We have assumed independent censoring in this paper. Another way in which patients may be treated 
differently following the occurrence of a clinically important event is to be withdrawn from a study. The 
occurrence of one event may increase the risk an investigator may withdraw the patient from the study 
and result in response-dependent censoring. If the events are independent conditional on the treatment 
covariate, this will not pose a problem but otherwise will lead to biased estimates of the baseline hazard 
functions and treatment effects. Use of inverse probability of censoring weights will help reduce this 
bias, and this is currently under investigation. 

Finally, we have focused on the frequency properties of estimators under a Cox regression models. 
Cox models are used routinely, but of course, the proportional hazards assumption may not be valid for 
either the marginal distributions or the composite endpoint. There is increasing interest in use of alter- 
native regression models for the analysis of survival data including accelerated failure time models and 
additive models. Exploration of the behavior of estimators from such models warrants study. 

Appendix: Derivation of the limiting value /? 

Suppose the proportional hazards assumption holds for the marginal distribution of each component 
event and a copula model is used to characterize the association. Under a working independence 
assumption [20], the limiting value for fi k is /3 k . Here, we consider a generic scalar covariate Z,, 
but note that the functions simplify with Z, a binary treatment indicator because Z[ = Z, in 



this case. Let S?\t) = m~ l Y!?=\ Y ik (t)X k0 (t)exp(fi k Zi)Zf, s%\t) = E{S ( k r) (t)), S r k {fi k ,t) = 
m- l YT= i Y ik {t)^ V {fi k Z i )Z r i , and s ( k r) (p k ,t) = E{S^{fi k ,t)), r = 0,1,2, where E(-) denotes 



expectation with respect to the true distribution. Let A(fi) — diag{j4jfc(/Jjt)}, k — 1, . . . , K, where the 
kth. diagonal element of A(fi) is 
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by the Theorem 4.2 of Andersen and Gill [35]. In the present setting, the true model is known, 
the required expectations can be obtained in closed form, and the integral can be evaluated using 
numerical methods. 
If 

M ik (t) = N ik (t) - J Y ik (t)X ok (t)e^ z 'dt 
is the martingale for events of type k; let 

Wikifa) = / \Zi - k m H dM ik (t), 

and w,G6) = (w n (Pi), • • . WiK^K))'- Then if we define B(fi) = E(yri(fi)\fi(fi)'), the asymptotic 
robust covariance matrix £ (/3) takes the form A(fi)~ l B{P)A{fi)~ l [20]. This can be used to obtain the 
limiting value through (12). 

The entries of B(fi) are obtained as follows. The (J, j ) element of B(fi) is 

E(wfj(0j,t)) = E{w tJ (pj),w t j(Pj)) 

h \\ sfgtj.t)) ' ) 

where (•, •) is the predictable covariation process and the last equality holds because of Fubini's theorem 
[36]. The (J, k) element of B(fi) is then 

E(Wij(pj)w ik (P k )) = E{Wij(fij), w ik (fi k )) 

= E i ( z '-fe)( z '-^) ( ^ fe, - dMi(,i)> ' 

Using the covariance function for correlated martingales of Prentice and Cai [37], (dMj (tj), dM k (t k )) 
can be obtained. In the case of bivariate data, {dM\(t\), dM 2 (t 2 )) is obtained simply as 

{dM 1 (t 1 ),dM 2 (t 2 )) = J r (dt 1 ,dt 2 \z i : ^dt^ + J r (ti,dt 2 \Z i ; fi)Ai(dfi|z<; Q)dtidt 2 
+ T{dh,t 2 \zf, Q)A 2 (dt 2 \zi; 0)dtidt 2 
+ J r (dt l ,dt 2 \Zi;Sl)A 1 (dt 1 )A 2 (dt 2 \zi] £l)dhdt 2 , 

where A k (dt k \zi;^) = dA k (t k \zi;^)/dt k ; F(dh,dt 2 \zi;Si) = d 2 T(t l ,t 2 \z i \^)/dt 1 dt 2 , 
T(dh,t 2 \Zi;&) = dJ r (ti,t 2 \z.i;Q).3t u and JY/t,df 2 |z;; ft) = dT(h, t 2 \Z { ; Q)/dt 2 . More specif- 
ically, if the joint survivor function F(ti, t 2 \zt; £2) is specified by the Clayton copula with mar- 
gins of two exponential distributions, then (dM/(f/), dM k (t k )) can be obtained in closed form and 
E(wij(fi j)Wik(fik)) can De obtained through numerical integration. Thus, we obtain the limiting value 
of robust variance, then the limiting weights can be calculated using c(/3) = T-itfyj'T- 1 ! and the 
limiting value /3 using equation (12). 
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